!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This file contains the simulated annealing algorithm
! used in the estimation. See the end of the file for
! attribution and details.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    MODULE simulated_anneal

    IMPLICIT NONE

    INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60)

    ! The following variables were in COMMON /raset1/
    REAL, SAVE    :: u(97), cc, cd, cm
    INTEGER, SAVE :: i97, j97



    CONTAINS

    SUBROUTINE sa(poe,nparam, x0, max, rt, eps0, ns, nt, neps, maxevl, lb, ub, step, iprint,  &
        iseed1, iseed2, temper, vm, xopt, fopt, nacc, nfcnev, nobds, ier)


    !use parameters
    !use globals

    !  Type all external variables.
    real (8), INTENT(IN)     :: lb(:), ub(:), step(:), eps0, rt  !chang all c to step
    real (8), INTENT(IN OUT) :: x0(:), temper, vm(:)
    real (8), INTENT(OUT)    :: xopt(:), fopt
    INTEGER, INTENT(IN)       :: nparam, ns, nt, neps, maxevl, iprint, iseed1, iseed2, poe
    INTEGER, INTENT(OUT)      :: nacc, nfcnev, nobds, ier
    LOGICAL, INTENT(IN)       :: max
    !real*8, intent(in):: x0(nparam),lb(nparam),ub(nparam),vm(nparam),temper,step(nparam),rt,eps0

    !
    !logical::max


    !  Type all internal variables.
    real (8) :: f, fp, panneal, pp, ratio, xp(nparam), fstar(neps)

    INTEGER   :: nup, ndown, nrej, nnew, lnobds, h, i, j, m, nacp(nparam),nmom=7, nnsim
    LOGICAL   :: quit





    !  Initialize the random number generator RANMAR.
    call rmarin(iseed1, iseed2)

    !  Set initial values.
    nacc = 0
    nobds = 0
    nfcnev = 0
    ier = 99

    DO i = 1, nparam
        xopt(i) = x0(i)  !initialize xoptimal to initial parameter.
        nacp(i) = 0
    END DO

    fstar = 1.0D+20
    OPEN(unit=200,file='ANNEALING.txt');
    OPEN(unit=201,file='ANNEALING1.txt');
    !  If the initial temperature is not positive, notify the user and
    !  return to the calling routine.
    IF (temper <= 0.0) THEN
        WRITE(200,'(/, "  THE INITIAL TEMPERATURE IS NOT POSITIVE. "/,  &
            &    "  reset the variable temper. "/)')
        ier = 3
        RETURN
    END IF
    !print *, 'b'
    !  If the initial value is out of bounds, notify the user and return
    !  to the calling routine.
    DO i = 1, nparam

        IF ((x0(i) > ub(i)) .OR. (x0(i) < lb(i))) THEN   !when parameter intial is out of boundary
            !print *, 'c'
            !pause
            call prt1()
            ier = 2
            RETURN
        END IF
    END DO

    !  Evaluate the function with input x0 and return value as F.
    fopt = 100000000     !give fopt a very big number


    call fcn(nparam,x0,f,poe)       ! give x0, call fcn,which solve for policy,simulate model and calculate loss function
    !print *, 'd'


    !  If the function is to be minimized, switch the sign of the function.
    !  Note that all intermediate and final output switches the sign back
    !  to eliminate any possible confusion for the user.
    IF(.NOT. max) f = -f
    nfcnev = nfcnev + 1
    fopt = f
    fstar(1) = f
    IF(iprint >= 1)  call prt2(max, nparam, x0, f)

    !  Start the main loop. Note that it terminates if (i) the algorithm
    !  succesfully optimizes the function or (ii) there are too many
    !  function evaluations (more than MAXEVL).

    nnsim=0
100 nup = 0
    nrej = 0
    nnew = 0
    ndown = 0
    lnobds = 0

    DO m = 1, nt
        DO j = 1, ns
            DO h = 1, nparam

                !  Generate XP, the trial value of x0. Note use of VM to choose XP.
                DO i = 1, nparam
                    IF (i == h) THEN
                        xp(i) = x0(i) + (ranmar()*2. - 1.) * vm(i)
                    ELSE
                        xp(i) = x0(i)
                    END IF

                    !  If XP is out of bounds, select a point in bounds for the trial.
                    IF((xp(i) < lb(i)) .OR. (xp(i) > ub(i))) THEN
                        !print *,xp
                        !   print *, 'hey1'
                        !pause
                        xp(i) = lb(i) + (ub(i) - lb(i))*ranmar()
                        !    print *,xp
                        !   print *, 'hey'
                        !pause
                        lnobds = lnobds + 1
                        nobds = nobds + 1
                        IF(iprint >= 3)  call prt3(max, nparam, xp, x0, f)
                        ! print *, 'e'
                        ! pause
                    END IF
                END DO


                !  Evaluate the function with the trial point XP and return as FP.
                !print *,xp

                !print *, fp
                !pause
                call fcn(nparam,xp,fp,poe)
                !   print *, 'f'
                !  pause
                IF(.NOT. max) fp = -fp
                nfcnev = nfcnev + 1

                IF(iprint >= 3)  call prt4(max, nparam, xp, x0, fp, f)

                !  If too many function evaluations occur, terminate the algorithm.
                IF(nfcnev >= maxevl) THEN
                    call prt5()
                    IF (.NOT. max) fopt = -fopt
                    ier = 1
                    RETURN
                END IF

                !  Accept the new point if the function value increases.
                IF(fp >= f) THEN
                    IF(iprint >= 3) THEN
                        WRITE(200,'("  POINT ACCEPTED")')

                    END IF
                    x0(1:nparam) = xp(1:nparam)
                    f = fp
                    nacc = nacc + 1
                    nacp(h) = nacp(h) + 1
                    nup = nup + 1

                    !  If greater than any other point, record as new optimum.
                    IF (fp > fopt) THEN
                        IF(iprint >= 3) THEN
                            WRITE(200,'("  NEW OPTIMUM")')
                            ! print *, 'new optimum'
                            !  call prt8(nparam, vm, xopt, x0)
                        END IF

                        xopt(1:nparam) = xp(1:nparam)
                        fopt = fp

                        print *, 'The bestval so far is', -fopt
                        print *,'The optimal parameters are',xopt
                        !           WRITE(201,'("The bestval so far is")')
                        WRITE(201, '("  The bestval so far is:  ",G25.18)') -fopt

                        !WRITE(201,*),xopt
                        call prtvec1(xopt, nparam, 'The optimal parameters are')


                        nnew = nnew + 1
                    END IF

                    !  If the point is lower, use the Metropolis criteria to decide on
                    !  acceptance or rejection.
                ELSE
                    panneal = exprep((fp - f)/temper)
                    pp = ranmar()
                    IF (pp < panneal) THEN
                        IF(iprint >= 3)  call prt6(max)
                        x0(1:nparam) = xp(1:nparam)
                        f = fp
                        nacc = nacc + 1
                        nacp(h) = nacp(h) + 1
                        ndown = ndown + 1
                    ELSE
                        nrej = nrej + 1
                        IF(iprint >= 3)  call prt7(max)
                    END IF
                END IF

            END DO
        END DO

        !  Adjust VM so that approximately half of all evaluations are accepted.
        DO i = 1, nparam
            ratio = DBLE(nacp(i)) /DBLE(ns)
            IF (ratio > .6) THEN
                vm(i) = vm(i)*(1. + step(i)*(ratio - .6)/.4)
            ELSE IF (ratio < .4) THEN
                vm(i) = vm(i)/(1. + step(i)*((.4 - ratio)/.4))
            END IF
            IF (vm(i) > (ub(i)-lb(i))) THEN
                vm(i) = ub(i) - lb(i)
            END IF
        END DO

        IF(iprint >= 2) THEN
            call prt8(nparam, vm, xopt, x0)
        END IF

        nacp(1:nparam) = 0

    END DO

    IF(iprint >= 1) THEN
        call prt9(max,nparam,temper,xopt,vm,fopt,nup,ndown,nrej,lnobds,nnew,&
            nfcnev,nmom)
    END IF

    !  Check termination criteria.
    quit = .false.
    fstar(1) = f
    IF ((fopt - fstar(1)) <= eps0) quit = .true.
    DO i = 1, neps
        IF (ABS(f - fstar(i)) > eps0) quit = .false.
    END DO

    !  Terminate SA if appropriate.
    IF (quit) THEN
        x0(1:nparam) = xopt(1:nparam)
        ier = 0
        IF (.NOT. max) fopt = -fopt
        IF(iprint >= 1)  call prt10()
        RETURN
    END IF

    !  If termination criteria is not met, prepare for another loop.
    temper = rt*temper
    DO i = neps, 2, -1
        fstar(i) = fstar(i-1)
    END DO
    f = fopt
    x0(1:nparam) = xopt(1:nparam)




    !  Loop again.
    !nnsim=nnsim+1
    !print *, 'this is interation', nnsim
    GO TO 100

    END SUBROUTINE sa


    FUNCTION exprep(rdum) RESULT(fn_val)
    !  This function replaces exp to avoid under- and overflows and is
    !  designed for IBM 370 type machines. It may be necessary to modify
    !  it for other machines. Note that the maximum and minimum values of
    !  EXPREP are such that they has no effect on the algorithm.

    real (8), INTENT(IN) :: rdum
    real (8)             :: fn_val

    IF (rdum > 174._dp) THEN
        fn_val = 3.69D+75
    ELSE IF (rdum < -180.+dp) THEN
        fn_val = 0.0_dp
    ELSE
        fn_val = EXP(rdum)
    END IF

    RETURN
    END FUNCTION exprep


    SUBROUTINE rmarin(ij, kl)
    !  This subroutine and the next function generate random numbers. See
    !  the comments for SA for more information. The only changes from the
    !  orginal code is that (1) the test to make sure that RMARIN runs first
    !  was taken out since SA assures that this is done (this test didn'temper
    !  compile under IBM's VS Fortran) and (2) typing ivec as integer was
    !  taken out since ivec isn'temper used. With these exceptions, all following
    !  lines are original.

    ! This is the initialization routine for the random number generator
    !     RANMAR()
    ! NOTE: The seed variables can have values between:    0 <= IJ <= 31328
    !                                                      0 <= KL <= 30081

    INTEGER, INTENT(IN) :: ij, kl

    INTEGER :: i, j, k, l, ii, jj, m
    REAL    :: s, temper

    IF( ij < 0  .OR.  ij > 31328  .OR. kl < 0  .OR.  kl > 30081 ) THEN
        WRITE(200, '(A)') ' The first random number seed must have a value ',  &
            'between 0 AND 31328'
        WRITE(200, '(A)') ' The second seed must have a value between 0 and 30081'
        STOP
    END IF
    i = MOD(ij/177, 177) + 2
    j = MOD(ij    , 177) + 2
    k = MOD(kl/169, 178) + 1
    l = MOD(kl,     169)
    DO ii = 1, 97
        s = 0.0
        temper = 0.5
        DO jj = 1, 24
            m = MOD(MOD(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = MOD(53*l+1, 169)
            IF (MOD(l*m, 64) >= 32) THEN
                s = s + temper
            END IF
            temper = 0.5 * temper
        END DO
        u(ii) = s
    END DO
    cc = 362436.0 / 16777216.0
    cd = 7654321.0 / 16777216.0
    cm = 16777213.0 /16777216.0
    i97 = 97
    j97 = 33
    RETURN
    END SUBROUTINE rmarin


    FUNCTION ranmar() RESULT(fn_val)
    REAL :: fn_val

    ! Local variable
    REAL :: uni

    uni = u(i97) - u(j97)
    IF( uni < 0.0 ) uni = uni + 1.0
    u(i97) = uni
    i97 = i97 - 1
    IF(i97 == 0) i97 = 97
    j97 = j97 - 1
    IF(j97 == 0) j97 = 97
    cc = cc - cd
    IF( cc < 0.0 ) cc = cc + cm
    uni = uni - cc
    IF( uni < 0.0 ) uni = uni + 1.0
    fn_val = uni

    RETURN
    END FUNCTION ranmar


    SUBROUTINE prt1()
    !  This subroutine prints intermediate output, as does PRT2 through
    !  PRT10. Note that if SA is minimizing the function, the sign of the
    !  function value and the directions (up/down) are reversed in all
    !  output to correspond with the actual function optimization. This
    !  correction is because SA was written to maximize functions and
    !  it minimizes by maximizing the negative a function.

    WRITE(200, '(/, "  THE STARTING VALUE (x0) IS OUTSIDE THE BOUNDS "/,  &
        &     "  (lb AND ub). execution terminated without any"/,  &
        &     "  optimization. respecify x0, ub OR lb so that  "/,  &
        &     "  lb(i) < x0(i) < ub(i), i = 1, nparam. "/)')

    RETURN
    END SUBROUTINE prt1


    SUBROUTINE prt2(max, nparam, x0, f)

    real (8), INTENT(IN) ::  x0(:), f
    INTEGER, INTENT(IN)   ::  nparam
    LOGICAL, INTENT(IN)   ::  max

    WRITE(200, '("  ")')
    call prtvec(x0,nparam,'INITIAL x0')
    IF (max) THEN
        WRITE(200, '("  INITIAL F: ",/, G25.18)') f
    ELSE
        WRITE(200, '("  INITIAL F: ",/, G25.18)') -f
    END IF

    RETURN
    END SUBROUTINE prt2


    SUBROUTINE prt3(max, nparam, xp, x0, f)

    real (8), INTENT(IN) ::  xp(:), x0(:), f
    INTEGER, INTENT(IN)   ::  nparam
    LOGICAL, INTENT(IN)   ::  max

    WRITE(200, '("  ")')
    call prtvec(x0, nparam, 'CURRENT x0')
    IF (max) THEN
        WRITE(200, '("  CURRENT F: ", G25.18)') f
    ELSE
        WRITE(200, '("  CURRENT F: ", G25.18)') -f
    END IF
    call prtvec(xp, nparam, 'TRIAL x0')
    WRITE(200, '("  POINT REJECTED SINCE OUT OF BOUNDS")')

    RETURN
    END SUBROUTINE prt3


    SUBROUTINE prt4(max, nparam, xp, x0, fp, f)

    real (8), INTENT(IN) ::  xp(:), x0(:), fp, f
    INTEGER, INTENT(IN)   ::  nparam
    LOGICAL, INTENT(IN)   ::  max

    WRITE(200,'("  ")')
    call prtvec(x0,nparam,'CURRENT x0')
    IF (max) THEN
        WRITE(200,'("  CURRENT F: ",G25.18)') f
        call prtvec(xp,nparam,'TRIAL x0')
        WRITE(200,'("  RESULTING F: ",G25.18)') fp
    ELSE
        WRITE(200,'("  CURRENT F: ",G25.18)') -f
        call prtvec(xp,nparam,'TRIAL x0')
        WRITE(200,'("  RESULTING F: ",G25.18)') -fp
    END IF

    RETURN
    END SUBROUTINE prt4


    SUBROUTINE prt5()

    WRITE(200, '(/, "  TOO MANY FUNCTION EVALUATIONS; CONSIDER "/,  &
        &     "  increasing maxevl OR eps0, OR decreasing "/,  &
        &     "  nt OR rt. these results are likely TO be "/, "  poor.",/)')

    RETURN
    END SUBROUTINE prt5


    SUBROUTINE prt6(max)
    LOGICAL, INTENT(IN) ::  max

    IF (max) THEN
        WRITE(200,'("  THOUGH LOWER, POINT ACCEPTED")')
    ELSE
        WRITE(200,'("  THOUGH HIGHER, POINT ACCEPTED")')
    END IF

    RETURN
    END SUBROUTINE prt6


    SUBROUTINE prt7(max)
    LOGICAL, INTENT(IN) :: max

    IF (max) THEN
        WRITE(200,'("  LOWER POINT REJECTED")')
    ELSE
        WRITE(200,'("  HIGHER POINT REJECTED")')
    END IF

    RETURN
    END SUBROUTINE prt7


    SUBROUTINE prt8(nparam, vm, xopt, x0)  ! ,krange,brange

    real (8), INTENT(IN) :: vm(:), xopt(:), x0(:) !,krange(:),brange(:)
    INTEGER, INTENT(IN)   :: nparam

    WRITE(200,'(/, " intermediate results after step length adjustment", /)')
    call prtvec(vm, nparam, 'NEW STEP LENGTH (VM)')
    call prtvec(xopt, nparam, 'CURRENT OPTIMAL x0')

    call prtvec(x0, nparam, 'CURRENT x0')
    WRITE(200,'(" ")')

    RETURN
    END SUBROUTINE prt8


    SUBROUTINE prt9(max, nparam, temper, xopt, vm, fopt, nup, ndown, nrej, lnobds, nnew, &
        nfcnev, nmom)

    real (8), INTENT(IN) :: xopt(:), vm(:), temper, fopt !,krange(:),brange(:)
    INTEGER, INTENT(IN)   :: nparam, nup, ndown, nrej, lnobds, nnew, nfcnev, nmom
    LOGICAL, INTENT(IN)   :: max

    ! Local variable
    INTEGER :: totmov

    totmov = nup + ndown + nrej

    WRITE(200,'(/," intermediate results before next temperature reduction",/)')
    WRITE(200,'("  CURRENT TEMPERATURE:            ",G12.5)') temper
    IF (max) THEN
        WRITE(200, '("  MAX FUNCTION VALUE SO FAR:  ",G25.18)') fopt
        WRITE(200, '("  TOTAL MOVES:                ",I8)') totmov
        WRITE(200, '("     UPHILL:                  ",I8)') nup
        WRITE(200, '("     ACCEPTED DOWNHILL:       ",I8)') ndown
        WRITE(200, '("     REJECTED DOWNHILL:       ",I8)') nrej
        WRITE(200, '("  OUT OF BOUNDS TRIALS:       ",I8)') lnobds
        WRITE(200, '("  NEW MAXIMA THIS TEMPERATURE:",I8)') nnew
        WRITE(200, '("  TOTAL NUMBER OF EVALUATIONS:",I8)') nfcnev
    ELSE
        WRITE(200, '("  MIN FUNCTION VALUE SO FAR:  ",G25.18)') -fopt
        WRITE(200, '("  TOTAL MOVES:                ",I8)') totmov
        WRITE(200, '("     DOWNHILL:                ",I8)')  nup
        WRITE(200, '("     ACCEPTED UPHILL:         ",I8)')  ndown
        WRITE(200, '("     REJECTED UPHILL:         ",I8)')  nrej
        WRITE(200, '("  TRIALS OUT OF BOUNDS:       ",I8)')  lnobds
        WRITE(200, '("  NEW MINIMA THIS TEMPERATURE:",I8)')  nnew
        WRITE(200, '("  TOTAL NUMBER OF EVALUATIONS:",I8)') nfcnev
    END IF
    call prtvec(xopt, nparam, 'CURRENT OPTIMAL x0')

    call prtvec(vm, nparam, 'STEP LENGTH (VM)')
    WRITE(200, '(" ")')

    RETURN
    END SUBROUTINE prt9


    SUBROUTINE prt10()

    WRITE(200, '(/, "  SA ACHIEVED TERMINATION CRITERIA. IER = 0. ",/)')

    RETURN
    END SUBROUTINE prt10


    SUBROUTINE prtvec(vector, ncols, name)
    !  This subroutine prints the double precision vector named VECTOR.
    !  Elements 1 thru NCOLS will be printed. NAME is a character variable
    !  that describes VECTOR. Note that if NAME is given in the call to
    !  PRTVEC, it must be enclosed in quotes. If there are more than 10
    !  elements in VECTOR, 10 elements will be printed on each line.

    INTEGER, INTENT(IN)           :: ncols
    real (8), INTENT(IN)         :: vector(ncols)
    CHARACTER (LEN=*), INTENT(IN) :: name

    INTEGER :: i, lines, ll

    WRITE(200,1001) NAME

    IF (ncols > 10) THEN
        lines = INT(ncols/10.)

        DO i = 1, lines
            ll = 10*(i - 1)
            WRITE(200,1000) vector(1+ll:10+ll)
        END DO

        WRITE(200,1000) vector(11+ll:ncols)
    ELSE
        WRITE(200,1000) vector(1:ncols)
    END IF

1000 FORMAT( 10(g12.5, ' '))
1001 FORMAT(/, 25(' '), a)

    RETURN
    END SUBROUTINE prtvec



    SUBROUTINE prtvec1(vector, ncols, name)
    !  This subroutine prints the double precision vector named VECTOR.
    !  Elements 1 thru NCOLS will be printed. NAME is a character variable
    !  that describes VECTOR. Note that if NAME is given in the call to
    !  PRTVEC, it must be enclosed in quotes. If there are more than 10
    !  elements in VECTOR, 10 elements will be printed on each line.

    INTEGER, INTENT(IN)           :: ncols
    real (8), INTENT(IN)         :: vector(ncols)
    CHARACTER (LEN=*), INTENT(IN) :: name

    INTEGER :: i, lines, ll

    WRITE(201,1001) NAME

    IF (ncols > 10) THEN
        lines = INT(ncols/10.)

        DO i = 1, lines
            ll = 10*(i - 1)
            WRITE(201,1000) vector(1+ll:10+ll)
        END DO

        WRITE(201,1000) vector(11+ll:ncols)
    ELSE
        WRITE(201,1000) vector(1:ncols)
    END IF

1000 FORMAT( 10(g12.5, ' '))
1001 FORMAT(/, 25(' '), a)

    RETURN
    END SUBROUTINE prtvec1

    SUBROUTINE prtvec11(momopt, nmom)
    real(8), intent(in) :: momopt(:)
    integer, intent(in) :: nmom

    integer :: jj

    do jj=1,nmom
        write(222, "(F20.10)") momopt(jj)
    enddo
    write(222, "( )")

    END SUBROUTINE


    END MODULE simulated_anneal


    !  Version: 3.2
    !  Date: 1/22/94.
    !  Differences compared to Version 2.0:
    !     1. If a trial is out of bounds, a point is randomly selected
    !        from LB(i) to UB(i). Unlike in version 2.0, this trial is
    !        evaluated and is counted in acceptances and rejections.
    !        All corresponding documentation was changed as well.
    !  Differences compared to Version 3.0:
    !     1. If VM(i) > (UB(i) - LB(i)), VM is set to UB(i) - LB(i).
    !        The idea is that if temper is high relative to LB & UB, most
    !        points will be accepted, causing VM to rise. But, in this
    !        situation, VM has little meaning; particularly if VM is
    !        larger than the acceptable region. Setting VM to this size
    !        still allows all parts of the allowable region to be selected.
    !  Differences compared to Version 3.1:
    !     1. Test made to see if the initial temperature is positive.
    !     2. WRITE statements prettied up.
    !     3. References to paper updated.

    !  Synopsis:
    !  This routine implements the continuous simulated annealing global
    !  optimization algorithm described in Corana et al.'s article "Minimizing
    !  Multimodal Functions of Continuous Variables with the "Simulated Annealing"
    !  Algorithm" in the September 1987 (vol. 13, no. 3, pp. 262-280) issue of
    !  the ACM Transactions on Mathematical Software.

    !  A very quick (perhaps too quick) overview of SA:
    !     SA tries to find the global optimum of an nparam dimensional function.
    !  It moves both up and downhill and as the optimization process
    !  proceeds, it focuses on the most promising area.
    !     To start, it randomly chooses a trial point within the step length
    !  VM (a vector of length nparam) of the user selected starting point. The
    !  function is evaluated at this trial point and its value is compared
    !  to its value at the initial point.
    !     In a maximization problem, all uphill moves are accepted and the
    !  algorithm continues from that trial point. Downhill moves may be
    !  accepted; the decision is made by the Metropolis criteria. It uses temper
    !  (temperature) and the size of the downhill move in a probabilistic
    !  manner. The smaller temper and the size of the downhill move are, the more
    !  likely that move will be accepted. If the trial is accepted, the
    !  algorithm moves on from that point. If it is rejected, another point
    !  is chosen instead for a trial evaluation.
    !     Each element of VM periodically adjusted so that half of all
    !  function evaluations in that direction are accepted.
    !     A fall in temper is imposed upon the system with the RT variable by
    !  temper(i+1) = RT*temper(i) where i is the ith iteration. Thus, as temper declines,
    !  downhill moves are less likely to be accepted and the percentage of
    !  rejections rise. Given the scheme for the selection for VM, VM falls.
    !  Thus, as temper declines, VM falls and SA focuses upon the most promising
    !  area for optimization.

    !  The importance of the parameter temper:
    !     The parameter temper is crucial in using SA successfully. It influences
    !  VM, the step length over which the algorithm searches for optima. For
    !  a small intial temper, the step length may be too small; thus not enough
    !  of the function might be evaluated to find the global optima. The user
    !  should carefully examine VM in the intermediate output (set IPRINT =
    !  1) to make sure that VM is appropriate. The relationship between the
    !  initial temperature and the resulting step length is function
    !  dependent.
    !     To determine the starting temperature that is consistent with
    !  optimizing a function, it is worthwhile to run a trial run first. Set
    !  RT = 1.5 and temper = 1.0. With RT > 1.0, the temperature increases and VM
    !  rises as well. Then select the temper that produces a large enough VM.

    !  For modifications to the algorithm and many details on its use,
    !  (particularly for econometric applications) see Goffe, Ferrier
    !  and Rogers, "Global Optimization of Statistical Functions with
    !  Simulated Annealing," Journal of Econometrics, vol. 60, no. 1/2,
    !  Jan./Feb. 1994, pp. 65-100.
    !  For more information, contact
    !              Bill Goffe
    !              Department of Economics and International Business
    !              University of Southern Mississippi
    !              Hattiesburg, MS  39506-5072
    !              (601) 266-4484 (office)
    !              (601) 266-4920 (fax)
    !              bgoffe@whale.st.usm.edu (Internet)

    !  As far as possible, the parameters here have the same name as in
    !  the description of the algorithm on pp. 266-8 of Corana et al.

    !  In this description, SP is single precision, DP is double precision,
    !  INT is integer, L is logical and (nparam) denotes an array of length nparam.
    !  Thus, DP(nparam) denotes a double precision array of length nparam.

    !  Input Parameters:
    !    Note: The suggested values generally come from Corana et al. To
    !          drastically reduce runtime, see Goffe et al., pp. 90-1 for
    !          suggestions on choosing the appropriate RT and NT.
    !    nparam - Number of variables in the function to be optimized. (INT)
    !    x0 - The starting values for the variables of the function to be
    !        optimized. (DP(nparam))
    !    MAX - Denotes whether the function should be maximized or minimized.
    !          A true value denotes maximization while a false value denotes
    !          minimization.  Intermediate output (see IPRINT) takes this into
    !          account. (L)
    !    RT - The temperature reduction factor.  The value suggested by
    !         Corana et al. is .85. See Goffe et al. for more advice. (DP)
    !    eps0 - Error tolerance for termination. If the final function
    !          values from the last neps temperatures differ from the
    !          corresponding value at the current temperature by less than
    !          eps0 and the final function value at the current temperature
    !          differs from the current optimal function value by less than
    !          eps0, execution terminates and IER = 0 is returned. (EP)
    !    NS - Number of cycles.  After NS*nparam function evaluations, each element of
    !         VM is adjusted so that approximately half of all function evaluations
    !         are accepted.  The suggested value is 20. (INT)
    !    NT - Number of iterations before temperature reduction. After
    !         NT*NS*nparam function evaluations, temperature (temper) is changed
    !         by the factor RT.  Value suggested by Corana et al. is
    !         MAX(100, 5*nparam).  See Goffe et al. for further advice. (INT)
    !    NEPS - Number of final function values used to decide upon termi-
    !           nation.  See eps0.  Suggested value is 4. (INT)
    !    MAXEVL - The maximum number of function evaluations.  If it is
    !             exceeded, IER = 1. (INT)
    !    LB - The lower bound for the allowable solution variables. (DP(nparam))
    !    UB - The upper bound for the allowable solution variables. (DP(nparam))
    !         If the algorithm chooses x0(I) .LT. LB(I) or x0(I) .GT. UB(I),
    !         I = 1, nparam, a point is from inside is randomly selected. This
    !         This focuses the algorithm on the region inside UB and LB.
    !         Unless the user wishes to concentrate the search to a particular
    !         region, UB and LB should be set to very large positive
    !         and negative values, respectively.  Note that the starting
    !         vector x0 should be inside this region.  Also note that LB and
    !         UB are fixed in position, while VM is centered on the last
    !         accepted trial set of variables that optimizes the function.
    !    step - Vector that controls the step length adjustment.  The suggested
    !        value for all elements is 2.0. (DP(nparam))
    !    IPRINT - controls printing inside SA. (INT)
    !             Values: 0 - Nothing printed.
    !                     1 - Function value for the starting value and
    !                         summary results before each temperature
    !                         reduction. This includes the optimal
    !                         function value found so far, the total
    !                         number of moves (broken up into uphill,
    !                         downhill, accepted and rejected), the
    !                         number of out of bounds trials, the
    !                         number of new optima found at this
    !                         temperature, the current optimal x0 and
    !                         the step length VM. Note that there are
    !                         nparam*NS*NT function evalutations before each
    !                         temperature reduction. Finally, notice is
    !                         is also given upon achieveing the termination
    !                         criteria.
    !                     2 - Each new step length (VM), the current optimal
    !                         x0 (XOPT) and the current trial x0 (x0). This
    !                         gives the user some idea about how far x0
    !                         strays from XOPT as well as how VM is adapting
    !                         to the function.
    !                     3 - Each function evaluation, its acceptance or
    !                         rejection and new optima. For many problems,
    !                         this option will likely require a small tree
    !                         if hard copy is used. This option is best
    !                         used to learn about the algorithm. A small
    !                         value for MAXEVL is thus recommended when
    !                         using IPRINT = 3.
    !             Suggested value: 1
    !             Note: For a given value of IPRINT, the lower valued
    !                   options (other than 0) are utilized.
    !    ISEED1 - The first seed for the random number generator RANMAR.
    !             0 <= ISEED1 <= 31328. (INT)
    !    ISEED2 - The second seed for the random number generator RANMAR.
    !             0 <= ISEED2 <= 30081. Different values for ISEED1
    !             and ISEED2 will lead to an entirely different sequence
    !             of trial points and decisions on downhill moves (when
    !             maximizing).  See Goffe et al. on how this can be used
    !             to test the results of SA. (INT)

    !  Input/Output Parameters:
    !    temper - On input, the initial temperature. See Goffe et al. for advice.
    !        On output, the final temperature. (DP)
    !    VM - The step length vector. On input it should encompass the region of
    !         interest given the starting value x0.  For point x0(I), the next
    !         trial point is selected is from x0(I) - VM(I)  to  x0(I) + VM(I).
    !         Since VM is adjusted so that about half of all points are accepted,
    !         the input value is not very important (i.e. is the value is off,
    !         SA adjusts VM to the correct value). (DP(nparam))

    !  Output Parameters:
    !    XOPT - The variables that optimize the function. (DP(nparam))
    !    FOPT - The optimal value of the function. (DP)
    !    NACC - The number of accepted function evaluations. (INT)
    !    NFCNEV - The total number of function evaluations. In a minor
    !             point, note that the first evaluation is not used in the
    !             core of the algorithm; it simply initializes the
    !             algorithm. (INT).
    !    NOBDS - The total number of trial function evaluations that
    !            would have been out of bounds of LB and UB. Note that
    !            a trial point is randomly selected between LB and UB. (INT)
    !    IER - The error return number. (INT)
    !          Values: 0 - Normal return; termination criteria achieved.
    !                  1 - Number of function evaluations (NFCNEV) is
    !                      greater than the maximum number (MAXEVL).
    !                  2 - The starting value (x0) is not inside the
    !                      bounds (LB and UB).
    !                  3 - The initial temperature is not positive.
    !                  99 - Should not be seen; only used internally.

    !  Work arrays that must be dimensioned in the calling routine:
    !       RWK1 (DP(NEPS))  (FSTAR in SA)
    !       RWK2 (DP(nparam))     (XP    "  " )
    !       IWK  (INT(nparam))    (NACP  "  " )
    !  nparam.B. In the Fortran 90 version, these are automatic arrays.

    !  Required Functions (included):
    !    EXPREP - Replaces the function EXP to avoid under- and overflows.
    !             It may have to be modified for non IBM-type main-
    !             frames. (DP)
    !    RMARIN - Initializes the random number generator RANMAR.
    !    RANMAR - The actual random number generator. Note that
    !             RMARIN must run first (SA does this). It produces uniform
    !             random numbers on [0,1]. These routines are from
    !             Usenet's comp.lang.fortran. For a reference, see
    !             "Toward a Universal Random Number Generator"
    !             by George Marsaglia and Arif Zaman, Florida State
    !             University Report: FSU-SCRI-87-50 (1987).
    !             It was later modified by F. James and published in
    !             "A Review of Pseudo-random Number Generators." For
    !             further information, contact stuart@ads.com. These
    !             routines are designed to be portable on any machine
    !             with a 24-bit or more mantissa. I have found it produces
    !             identical results on a IBM 3081 and a Cray Y-MP.

    !  Required Subroutines (included):
    !    PRTVEC - Prints vectors.
    !    PRT1 ... PRT10 - Prints intermediate output.
    !    FCN - Function to be optimized. The form is
    !            SUBROUTINE FCN(nparam, x0, F)
    !            IMPLICIT NONE
    !            INTEGER, PARAMETER     :: dp = SELECTED_REAL_KIND(14, 60)
    !            INTEGER, INTENT(IN)    :: nparam
    !            REAL (dp), INTENT(IN)  :: x0(nparam)
    !            REAL (dp), INTENT(OUT) :: F
    !            ...
    !            function code with F = F(x0)
    !            ...
    !            RETURN
    !            END
    !  Note: This is the same form used in the multivariable
    !        minimization algorithms in the IMSL edition 10 library.

    !  Machine Specific Features:
    !    1. EXPREP may have to be modified if used on non-IBM type main-
    !       frames. Watch for under- and overflows in EXPREP.
    !    2. Some FORMAT statements use G25.18; this may be excessive for
    !       some machines.
    !    3. RMARIN and RANMAR are designed to be protable; they should not
    !       cause any problems.

    ! We need this to access the file names


    ! ABSTRACT:
    !   Simulated annealing is a global optimization method that distinguishes
    !   between different local optima.  Starting from an initial point, the
    !   algorithm takes a step and the function is evaluated. When minimizing a
    !   function, any downhill step is accepted and the process repeats from this
    !   new point. An uphill step may be accepted. Thus, it can escape from local
    !   optima. This uphill decision is made by the Metropolis criteria. As the
    !   optimization process proceeds, the length of the steps decline and the
    !   algorithm closes in on the global optimum. Since the algorithm makes very
    !   few assumptions regarding the function to be optimized, it is quite
    !   robust with respect to non-quadratic surfaces. The degree of robustness
    !   can be adjusted by the user. In fact, simulated annealing can be used as
    !   a local optimizer for difficult functions.

    !   This implementation of simulated annealing was used in "Global Optimizatio
    !   of Statistical Functions with Simulated Annealing," Goffe, Ferrier and
    !   Rogers, Journal of Econometrics, vol. 60, no. 1/2, Jan./Feb. 1994, pp.
    !   65-100. Briefly, we found it competitive, if not superior, to multiple
    !   restarts of conventional optimization routines for difficult optimization
    !   problems.

    !   For more information on this routine, contact its author:
    !   Bill Goffe, bgoffe@whale.st.usm.edu

    ! This version in Fortran 90 has been prepared by Alan Miller.
    ! It is compatible with Lahey's ELF90 compiler.
    ! nparam.B. The 3 last arguments have been removed from subroutine sa.   these
    !      were work arrays and are now internal to the routine.
    ! e-mail:  amiller @ bigpond.net.au
    ! URL   :  http://users.bigpond.net.au/amiller

    ! Latest revision of Fortran 90 version - 3 August 1997
